AQUA and TERRA
Two satellites
So far, I have only been looking at data from the AQUA satellite. MAIAC AOT data comes from two different satellites: AQUA and TERRA. The quality of the data should be pretty much equivalent, but the time that data is collected is staggered by a few hours so comparing data from AQUA and TERRA should give us an idea of the rate at which conditions change. This will also help us determine how representative a measurement at a particular time might be.
oct08A <- maiac_loadRaster("h01v04", 20171008, 2105, converterPath = "../executables/h4toncff_nc4") %>%
projectRaster(crs = CRS("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0")) %>%
crop(raster::extent(-125.52, -117.03, 36.52, 41.78))## /Users/helen/Data/maiac/MAIACAAOT.h01v04.20172812105.nc already exists. Skipping conversionoct08T <- maiac_loadRaster("h01v04", 20171008, 1930, product = "MAIACTAOT", converterPath = "../executables/h4toncff_nc4") %>%
projectRaster(crs = CRS("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0")) %>%
crop(raster::extent(-125.52, -117.03, 36.52, 41.78))## /Users/helen/Data/maiac/MAIACTAOT.h01v04.20172811930.nc already exists. Skipping conversionoct09A <- maiac_loadRaster("h01v04", 20171009, 2150, converterPath = "../executables/h4toncff_nc4") %>%
projectRaster(crs = CRS("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0")) %>%
crop(raster::extent(-125.52, -117.03, 36.52, 41.78))## /Users/helen/Data/maiac/MAIACAAOT.h01v04.20172822150.nc already exists. Skipping conversionoct09T <- maiac_loadRaster("h01v04", 20171009, product = "MAIACTAOT", converterPath = "../executables/h4toncff_nc4") %>%
projectRaster(crs = CRS("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0")) %>%
crop(raster::extent(-125.52, -117.03, 36.52, 41.78))## /Users/helen/Data/maiac/MAIACTAOT.h01v04.20172821835.nc already exists. Skipping conversionoct10A <- maiac_loadRaster("h01v04", 20171010, converterPath = "../executables/h4toncff_nc4") %>%
projectRaster(crs = CRS("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0")) %>%
crop(raster::extent(-125.52, -117.03, 36.52, 41.78))## /Users/helen/Data/maiac/MAIACAAOT.h01v04.20172832055.nc already exists. Skipping conversionoct10T <- maiac_loadRaster("h01v04", 20171010, 1915, product = "MAIACTAOT", converterPath = "../executables/h4toncff_nc4") %>%
projectRaster(crs = CRS("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0")) %>%
crop(raster::extent(-125.52, -117.03, 36.52, 41.78))## /Users/helen/Data/maiac/MAIACTAOT.h01v04.20172831915.nc already exists. Skipping conversionoct11A <- maiac_loadRaster("h01v04", 20171011, 2135, converterPath = "../executables/h4toncff_nc4") %>%
projectRaster(crs = CRS("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0")) %>%
crop(raster::extent(-125.52, -117.03, 36.52, 41.78))## /Users/helen/Data/maiac/MAIACAAOT.h01v04.20172842135.nc already exists. Skipping conversionoct11T <- maiac_loadRaster("h01v04", 20171011, 1820, product = "MAIACTAOT", converterPath = "../executables/h4toncff_nc4") %>%
projectRaster(crs = CRS("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0")) %>%
crop(raster::extent(-125.52, -117.03, 36.52, 41.78))## /Users/helen/Data/maiac/MAIACTAOT.h01v04.20172841820.nc already exists. Skipping conversionoct12A <- maiac_loadRaster("h01v04", 20171012, 2040, converterPath = "../executables/h4toncff_nc4") %>%
projectRaster(crs = CRS("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0")) %>%
crop(raster::extent(-125.52, -117.03, 36.52, 41.78))## /Users/helen/Data/maiac/MAIACAAOT.h01v04.20172852040.nc already exists. Skipping conversionoct12T <- maiac_loadRaster("h01v04", 20171012, product = "MAIACTAOT", converterPath = "../executables/h4toncff_nc4") %>%
projectRaster(crs = CRS("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0")) %>%
crop(raster::extent(-125.52, -117.03, 36.52, 41.78))## /Users/helen/Data/maiac/MAIACTAOT.h01v04.20172851905.nc already exists. Skipping conversionoct13A <- maiac_loadRaster("h01v04", 20171013, converterPath = "../executables/h4toncff_nc4") %>%
projectRaster(crs = CRS("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0")) %>%
crop(raster::extent(-125.52, -117.03, 36.52, 41.78))## /Users/helen/Data/maiac/MAIACAAOT.h01v04.20172862125.nc already exists. Skipping conversionoct13T <- maiac_loadRaster("h01v04", 20171013, 1945, product = "MAIACTAOT", converterPath = "../executables/h4toncff_nc4") %>%
projectRaster(crs = CRS("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0")) %>%
crop(raster::extent(-125.52, -117.03, 36.52, 41.78))## /Users/helen/Data/maiac/MAIACTAOT.h01v04.20172861945.nc already exists. Skipping conversionSide-by-side
colors <- colorRampPalette(c("white", "orange", "red"))(9)
breaks <- c(0, 0.03, 0.1, 0.2, 0.3, 0.5, 0.8, 10)
rasterList <- list(oct08A, oct08T, oct09A, oct09T, oct10A, oct10T, oct11A, oct11T, oct12A, oct12T, oct13A, oct13T)
par(mfcol = c(1,2))
for ( i in 1:6 ) {
AQUA <- rasterList[[i*2 - 1]]
TERRA <- rasterList[[i*2]]
day <- stringr::str_pad(i+7, 2, "left", "0")
plot(AQUA, breaks = breaks, col = colors, colNA = "gray90", alpha = .8, main = paste0("Oct ", day), legend = FALSE)
plot(USCensusStates, add = TRUE)
plot(TERRA, breaks = breaks, col = colors, colNA = "gray90", alpha = .8, main = paste0("Oct ", day))
plot(USCensusStates, add = TRUE)
}
Difference
How does it change from TERRA (which is generally earlier) to AQUA?
colors <- colorRampPalette(c("blue", "white", "red"))(13)
pow <- c(-2, -1.5, -1, -.5, 0, .5, 1)
breaks <- c(-10^rev(pow), 10^pow)
par(mar = c(3.1, 3.1, 4.1, 4.1))
for ( i in 1:6 ) {
AQUA <- rasterList[[i*2 - 1]]
TERRA <- rasterList[[i*2]]
day <- stringr::str_pad(i+7, 2, "left", "0")
plot(AQUA-TERRA, col = colors, breaks = breaks, colNA = "gray90", alpha = .8, main = paste0("Oct ", day), legend = FALSE)
plot(USCensusStates, add = TRUE)
usr <- par("usr")
legendbrks <- seq(usr[3], usr[4], length = length(breaks))
lx <- usr[2] + .5
rx <- usr[2] + .75
rect(lx,
head(legendbrks, -1),
rx,
tail(legendbrks, -1),
col = colors,
xpd = NA)
axis(4, at = legendbrks, labels = c(paste0("-10^", rev(pow)), paste0("10^", pow)), pos = rx, las = 2)
}
Would it make sense to combine them somehow to fill in missing values? Here are some possible methods:
colors <- colorRampPalette(c("white", "orange", "red"))(9)
breaks <- c(0, 0.03, 0.1, 0.2, 0.3, 0.5, 0.8, 3)
AQUA <- oct10A
TERRA <- oct10T
plot(max(AQUA, TERRA, na.rm = TRUE), breaks = breaks, col = colors, colNA = "gray90", main = "Oct 10: max values")
plot(USCensusStates, add = TRUE)
plot(mean(AQUA, TERRA, na.rm = TRUE), breaks = breaks, col = colors, colNA = "gray90", main = "Oct 10: mean values")
plot(USCensusStates, add = TRUE)
plot(AQUA, breaks = breaks, col = colors, colNA = "gray90", main = "Oct 10: TERRA over AQUA")
plot(TERRA, breaks = breaks, col = colors, add = TRUE)
plot(USCensusStates, add =TRUE)
plot(TERRA, breaks = breaks, col = colors, colNA = "gray90", main = "Oct 10: AQUA over TERRA")
plot(AQUA, breaks = breaks, col = colors, add = TRUE)
plot(USCensusStates, add = TRUE)